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Abstract 

Existing  models  for  predicting  the  penetration  depth  of  munitions  and  ex¬ 
plosives  of  concern  are  inaccurate  and  insufficient  from  a  user  (range 
manager,  U.S.  Army  Corps  of  Engineers  project  manager,  or  environmen¬ 
tal  consultant)  operability  perspective  for  current  needs.  We  attribute  poor 
model  performance  to  ( 1)  a  heavy  dependence  on  empirically  derived  pa- 
rameterizations  poorly  linked  to  the  physical  properties  of  the  target  mate¬ 
rial  or  (2)  physics- based  models  that  inadeguately  capture  the  salient  me¬ 
chanical  processes,  especially  in  the  first  meter  of  penetration.  Conse- 
guently,  we  have  developed  a  micromechanical-based  model  using  a  hy¬ 
brid  discrete  element  model  (DEM)  /  finite  element  model  (FEM)  ap¬ 
proach  capable  of  a  detailed  treatment  of  near- surface  soil  properties.  To 
examine  the  effects  of  varying  levels  of  moisture  on  the  dynamic  behavior 
of  a  soil,  we  fabricated  a  small-scale  biaxial  shear  test  to  inform  the  devel¬ 
opment  and  calibration  of  the  DEM  contact  model.  We  conducted  projec¬ 
tile-drop  tests  into  sand  with  a  scale  version  of  a  57  mm  projectile  and 
measured  projechle  penetration  to  compare  with  model  results. 


DISCLAIMER:  The  contents  of  this  report  are  not  to  be  used  for  advertising,  publication,  or  promotional  purposes.  Ci¬ 
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All  product  names  and  trademarks  cited  are  the  property  of  their  respective  owners.  The  findings  of  this  report  are  not  to 
be  construed  as  an  official  Department  of  the  Army  position  unless  so  designated  by  other  authorized  documents. 
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Executive  Summary 

Objectives 

Existing  models  for  predicting  penetration  depth  of  munitions  and  explo¬ 
sives  of  concern  (MEC)  are  inaccurate  and  insufficient  from  a  user  (range 
manager,  U.S.  Army  Corps  of  Engineers  (USACE)  project  manager,  or  en¬ 
vironmental  consultant)  operability  perspective  for  current  needs.  We  at¬ 
tribute  poor  model  performance  to  ( 1)  a  heavy  dependence  on  empirically 
derived  parameterizations  poorly  linked  to  the  physical  properties  of  the 
target  material  or  (2)  physics- based  models  that  inadeguately  capture  the 
salient  mechanical  processes,  especially  in  the  first  meter  of  penetration. 
To  address  these  shortcomings  our  objective  is  to  develop  improved  con¬ 
stitutive  behavior  models  by  using  a  micromechanical  approach  that  ex¬ 
plicitly  accounts  for  material  properties  such  as  soil  moisture  content, 
grain  size,  shape,  and  density. 

Technical  Approach 

Our  technical  approach  involves  the  development  of  a  micromechanical- 
based  model  using  a  hybrid  discrete  element  model  (DEM)  /  finite  ele¬ 
ment  model  (FEM)  capable  of  a  detailed  treatment  of  near- surface  soil 
properties.  DEM  model  particle  configurations  were  generated  from  three- 
dimensional  microCr  imaging  of  a  soil  sample.  To  examine  the  effects  of 
varying  levels  of  moisture  on  the  dynamic  behavior  of  a  soil,  we  fabricated 
a  small-scale  triaxial  shear  test  to  inform  the  DEM  contact  model  develop¬ 
ment  and  calibration.  Projectile- drop  tests  into  sand  were  conducted  with 
a  scale  version  of  a  57  mm  projectile  where  we  measured  projectile  pene¬ 
tration  for  comparison  with  model  results. 

Results 

An  improved  constitutive  model  framework  has  been  developed  that  im¬ 
proves  proj  ectile  penetration  by  relying  only  on  the  parent  material  char¬ 
acteristics  (elastic  modulus,  Poisson's  ratio),  grain  geometry,  friction  coef¬ 
ficient,  and  the  volume  of  pore  water.  Of  these  parameters,  the  friction  co¬ 
efficient  is  the  least  known  for  a  particular  soil  type.  The  triaxial  shear  test 
experiments,  which  characterize  the  dynamic  behavior  and  transfer  of  en¬ 
ergy  through  the  soil  fabric,  serve  as  a  good  source  of  data  to  calibrate  the 
friction  coefficient.  The  preliminary  numerical  model  results  compare  well 
in  a  gualitative  sense  with  the  drop- test  measurements. 


ERDC/CRREL  TR-17-12 


viii 


Benefits 

Application  of  these  improved  constitutive  soil  behavior  models  allows  for 
probability  estimates  of  specific  munition  types  to  the  depth  of  interest 
that  will  lead  to  accurate  time  and  cost  projections  of  MEC  cleanup  during 
project  planning.  This  will  provide  the  MEC  recovery  team  with  tighter 
bounds  on  the  uncertainty  associated  with  apparent  MEC  at  depth  identi¬ 
fied  during  geophysical  surveys  and  will  limit  the  over- excavation  of  ap¬ 
parent  MEC  at  depth,  thus  reducing  the  Department  of  Defense  cleanup 
costs  for  the  Formerly  Used  Defense  Sites  and  Base  Realignment  and  Clo¬ 
sure  sites  under  the  Military  Munition  Response  Program. 
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1  Introduction 

1.1  Background 

Over  4900  munitions  testing  sites  in  the  United  States  will  require  a  com¬ 
bined  total  of  $13  billion  of  remediation  efforts  over  the  following  decades 
to  remove  munitions  and  explosives  of  concern  (MEC).  The  munitions- 
range  managers,  U.S.  Army  Corps  of  Engineers  (USACE)  project  manag¬ 
ers,  and  environmental  consultants  tasked  with  these  remediation  efforts 
have  several  tools  at  their  disposal  to  detect  or  estimate  the  depth  of  MEC, 
including  geophysical  sensors  and  physics- based  numerical  models.  Alt¬ 
hough  these  tools  have  reduced  remediation  costs  significantly,  their  cur¬ 
rent  depth- prediction  capabilities  are  lacking.  Physics- based  models  are 
unable  to  capture  the  important  mechanical  process  of  penetration,  espe¬ 
cially  within  1  m  of  the  ground's  surface,  and  rely  heavily  on  empirical  pa- 
rameterizations  that  poorly  describe  the  target  materials'  physical  proper¬ 
ties.  The  uncertainty  associated  with  these  methods  can  lead  to  inflated 
time  and  cost  estimates  for  MEC  cleanup. 

1.2  Objectives 

This  study  was  funded  by  the  Strategic  Environmental  Research  and  De¬ 
velopment  Program  (SERDP)  and  was  conducted  in  support  of  the  FY17 
SERDP  SEED  call  for  proposal.  The  goal  of  this  study  was  to  address  these 
shortcomings  by  developing  improved  constitutive  behavior  models  from  a 
micromechanical  approach  that  explicitly  accounts  for  material  properties 
such  as  soil  moisture  content,  grain  size,  shape,  and  density.  Applying 
these  improved  constitutive  soil  behavior  models  allows  for  probability  es¬ 
timates  of  specific  munition  types  to  the  depth  of  interest  and  will  lead  to 
accurate  time  and  cost  projections  of  MEC  cleanup  during  project  plan¬ 
ning.  This  will  provide  the  MEC  recovery  team  with  tighter  bounds  on  the 
uncertainty  associated  with  apparent  MEC  at  depth  identified  during  geo¬ 
physical  surveys  and  will  limit  the  overexcavation  of  apparent  MEC  at 
depth,  thus  reducing  the  Department  of  Defense  cleanup  costs  for  the  For¬ 
merly  Used  Defense  Sites  and  Base  Realignment  and  Closure  sites  under 
the  Military  Munition  Response  Program. 
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1.3  Approach 

Our  technical  approach  involved  developing  a  micromechanical-based 
model  using  a  hybrid  discrete  element  model  (DEM)  /  finite  element 
model  (FEM)  capable  of  a  detailed  treatment  of  near- surface  soil  proper¬ 
ties.  DEM  model  particle  configurations  were  generated  from  three-di¬ 
mensional  microCr  imaging  of  a  soil  sample.  To  examine  the  effects  of 
varying  levels  of  moisture  on  the  dynamic  behavior  of  a  soil,  we  fabricated 
a  small-scale  triaxial  shear  test  to  inform  the  DEM  contact  model  develop¬ 
ment  and  calibration.  Projectile- drop  tests  into  sand  were  conducted  with 
a  scale  version  of  a  57  mm  projectile  where  we  measured  projectile  pene¬ 
tration  for  comparison  with  model  results. 
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2  Fundamental  Theory 

2.1  Penetration  dynamics 

The  dynamics  of  a  proj  ectile  penetrating  a  target  medium  is  a  function  of 
the  projectile  geometry,  projectile  impact  energy,  target  material  proper¬ 
ties,  and  the  contact  mechanics  between  the  projectile  and  the  target  mate¬ 
rial.  Previous  studies  have  primarily  explored  the  penetration- dynamics 
problem,  ultimately  to  estimate  penetration  trajectory  and  depth,  as  a  cav¬ 
ity-expansion  problem  where  the  projectile  pushes  the  target  material  out 
of  the  way  to  create  the  cavity  that  the  proj  ectile  then  occupies.  This  strain 
due  to  cavity  formation  induces  stresses  on  the  cavity  boundary  (i.e.,  the 
projectile-  soil  interface)  that  act  normal  (cavity  pressure)  and  tangentially 
(friction  between  projectile  and  target  material)  to  the  projectile  surface, 
decelerating  the  body  as  it  moves  through  the  target  medium  (Bless  et  al. 
2012;  Borget  al.  2012).  This  deceleration  can  be  described  generally  as 

dv  m) 

dt  ^  ' 

where  a,  /),  and  y  are  empirically  determined  constants,  as  with  Young's 
eguations  (Young  1967, 1997),  or  are  obtained  using  a  physics- based  ap¬ 
proach  that  uses  cavity  expansion  theory  (e.g.,  the  work  from  Forrestal  et 
al.  1992  and  more  recently  Shi  et  al.  2014). 

There  are  two  commonly  used  simplifications  of  the  general  penetration 
eguation:  (1)  the  Robins- Euler  simplification  that  assumes  the  projectile 
deceleration  is  velocity  independent,  therefore  a  =  0  and /)  =0,  and  (2)  the 
Poncelet  eguation  where  /)  =0  and  a  and  y  are  nonzero  positive  constants. 
The  Robins- Euler  eguation  applies  to  low  velocity  impacts  where  the 
boundary  location  between  plastic  and  elastic  deformation  regions  re¬ 
mains  nearly  constant  for  the  velocity  range  of  interest  because  the  stress 
front  propagates  much  more  guickly  than  the  projectile.  As  the  projectile 
impact  velocity  increases,  the  cavity  pressure  begins  to  exhibit  velocity  de¬ 
pendence;  and  the  deceleration  function  takes  the  form  that  Poncelet  pos¬ 
its.  This  form  of  the  general  penetration  eguation  is  the  one  that  is  typi¬ 
cally  used  to  derive  the  empirical  and  physics- based  penetration- depth 
eguations. 
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Several  researchers  have  attempted  to  empirically  correlate  the  cavity 
pressure  to  the  proj  ectile  impact  velocity  and  geometry  based  on  penetra¬ 
tion-depth  experiments  (Young  1967, 1997;  Bernard  1977, 1978;  Bernard 
and  Creighton  1978, 1979;  Allen  et  al.  1957).  These  relationships  are  not 
physics  based  and  reguire  well- characterized  experiments  to  characterize 
the  soil's  penetrability  index,  or  S- number.  The  generality  of  these  empiri¬ 
cally  based  correlations  is  limited  or  nonexistent  and  would  reguire  a  vast 
database  of  experiments  to  account  for  variations  in  soil  type,  stratigraphy, 
and  moisture  content.  In  addition,  a  statistical  picture  of  the  penetration 
depth  would  be  impossible  with  this  type  of  model  because  of  this  loose 
correlation  between  the  soil  properties  and  the  parametrization. 

2.2  Physics-based  penetration  models 

In  contrast  to  empirical- correlation- based  penetration  models,  physics- 
based  cavity  expansion  models  attempt  to  determine  the  cavity  pressure, 
or  penetration  resistance,  according  to  behavior  governed  by  soil  mechan¬ 
ics  models,  such  as  the  Mohr- Coulomb  theory.  A  fundamental  component 
of  cavity- expansion  theory  is  the  determination  of  the  propagation  velocity 
of  the  boundary  between  the  plastically  and  elastically  deforming  regions, 
typically  using  the  yield  surface  governed  by  the  Mohr- Coulomb  or 
Drucker-Prager  failure  criterion  and  assuming  that  cavity  expansion  is 
spherically  symmetric.  Studies  using  the  spherical  cavity  expansion  ap¬ 
proach  have  accurately  predicted  penetration  depth  in  porous  rock;  con¬ 
crete;  and  dry,  coarse  sands  (Forrestal  et  al.  1992;  Shi  et  al.  2014;  Forrestal 
1986) .  However,  these  theories  are  limited  to  ogive  or  conical- nose  config¬ 
urations  due  to  the  spherical- symmetry  assumption.  The  DEM  approach 
does  not  make  any  of  these  simplifying  assumptions  because  the  partide- 
to-particle  interaction  is  explidtly  described,  and  bulk  properties  such  as 
plastic  strain  are  determined  a  posteriori.  The  numerical  test  bed  allows  us 
to  interrogate  the  soil  sample  in  ways  that  are  extremely  difficult  to  repli¬ 
cate  in  physical  experiments  (e.g.,  measuring  the  strain  field  as  the  projec¬ 
tile  moves  through  the  target  medium) . 

2.2.1  Near-surface  soil  properties 

Soil  constitutive  modeling  is  based  primarily  on  Hooke's  law  of  linear  elas- 
tidty  and  Coulomb's  law  of  perfed  plastidty  at  small  strains  (e  <  10  ) .  Re¬ 
search  has  shown  that  soil  behavior  is  complex  and  nonlinear  and  exhibits 
anisotropic  time- dependent  behavior.  This  is  particularly  true  in  the  multi¬ 
phase  continuum  mechanics  of  unsaturated  soils  in  the  near  surface,  the 
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upper  1  m  of  soil.  The  soil  properties  and  skeletal  structure  of  the  near- sur¬ 
face  material  will  have  the  greatest  influence  on  penetration  depth  because 
these  physical  properties  significantly  affect  the  anisotropic  confinement 
of  the  soil  material.  Laboratory  investigations  (Hardin  and  Richart  1963; 
Richart  et  al.  1970;  Hardin  and  Dmevich  1972;  Iwasaki  et  al.  1974;  Tat- 
suoka  et  al  1978;  Seed  et  al.  1986)  in  stress  conditions  greater  than  the  up¬ 
per  1  m  of  overburden  have  shown  that  vertical  confining  stress  has  the 
most  significant  impact  on  soil  state  and  behavior.  Richart  et  al.  (1970) 
state  that  the  effects  of  effective  confining  pressure  acting  at  each  point  in 
a  soil  mass  depends  on  the  effective  overburden  pressures,  a0,  and  that 
compressional  and  shear  wave  velocities  are  dependent  on  (<r0)0-25.  There¬ 
fore,  any  effect  that  causes  a  change  in  a0  also  causes  a  change  in  wave- 
propagation  velocity,  skeletal  structure,  and  ultimately  soil  behavior.  How¬ 
ever,  the  validity  of  these  equations  and  soil  behavior  at  the  near- surface 
boundary  (i.e.,  the  upper  1  m  of  soil)  is  not  well  understood,  primarily  due 
to  the  limited  vertical  confining  stress,  the  high  degree  of  anisotropy,  and 
the  inapplicability  of  Mohr- Coulomb  soil  mechanics,  which  is  often  ig¬ 
nored  or  overgeneralized. 

2.2.2  Friction  coefficient 

As  described  above,  the  problem  of  projectile  penetration  in  soil  involves 
significant  target  material  deformation  and  particle  motion  in  the  near 
field  and  relatively  small  deformations  (nominally  elastic)  in  the  far  field. 
Both  of  these  regions  of  behavior  must  be  adequately  described  for  accu¬ 
rate  penetration- depth  estimation.  The  grain- to- grain  friction  coefficient, 
/r,  is  a  key  physical  quantity  that  governs  the  point  at  which  two  grains  will 
slide  relative  to  each  other  and  thus  establishes  the  stress  conditions  under 
which  extensive  straining  is  initiated.  Because  this  condition  essentially 
defines  the  boundary  between  near-  and  far- field  behaviors  in  the  soil,  n  is 
a  critically  important  quantity  in  the  present  effort.  Additionally,  the  fric¬ 
tion  between  the  soil  particles  and  the  metal  skin  of  the  penetrator  is  im¬ 
portant  to  the  analysis. 

Recent  work  (Cole  20 15)  has  demonstrated  with  experiments  that  certain 
types  of  naturally  occurring  grains  (generally  ones  with  smooth  surfaces, 
such  as  river  sand  or  other  weathered  materials)  exhibit  a  n  that  actually 
decreases  with  increasing  normal  force.  In  some  cases,  this  effect  can  be 
quite  pronounced,  decreasing  from  approximately  0.8  to  0.2  as  the  normal 
force  increases  from  1  to  20  N.  This  has  the  effect  of  making  large-scale 
(e.g.,  permanent)  strain  much  easier  at  high  background  stress  levels  and 
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is  thus  anticipated  to  be  an  important  effect  to  include  in  the  penetration 
model. 
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3  Materials  and  Methods 

3.1  Discrete  element  method 

Our  approach  determines  mechanical  properties  of  soil  in  terms  of 
measureable  physical  properties  (particle  size,  shape,  density,  and  mois¬ 
ture  conditions)  using  a  DEM  model  where  the  particle- scale  mechanics  of 
a  soil  are  explicitly  modeled.  This  physics- based  method  provides  the 
means  to  derive  a  soil's  bulk  macroscopic  properties  directly  from  a  micro¬ 
mechanical  description  of  the  target  material.  Modeling  the  penetration 
process  by  using  the  DEM  involves  calculating  the  contact  forces  and  re¬ 
sulting  motion  for  each  individual  soil  grain  (for  a  grain  diameter  of  1  mm, 
a  meter  cubed  volume  would  reguire  approximately  12  million  particles), 
which  is  very  computationally  burdensome,  to  generate  robust  statistics 
that  reguire  a  large  number  of  model  realizations.  This  report  describes 
the  preliminary  steps  to  use  a  DEM  method  to  improve  penetration  mod¬ 
els  through  the  improvement  of  constitutive  models  used  for  soils  at  the 
near  surface. 

We  use  a  mechanistic  approach  to  determine  the  constitutive  behavior  of 
the  target  material  by  explicitly  describing  the  interaction  between  grains 
composing  the  target  material  and  the  projectile.  Our  approach  takes  ad¬ 
vantage  of  the  extensive  work  conducted  in  recent  years  by  the  Engineer 
Research  and  Development  Center  (ERDC),  Cold  Regions  Research  and 
Engineering  Laboratory  (CRREL),  on  granular  media  mechanics  and  mod¬ 
eling.  This  approach  considers  near-  and  far- field  behavior  of  the  soil.  The 
near  field  is  defined  as  the  region  of  soil  subjected  to  impact  loading  and 
large  strains.  This  region  reguires  certain  types  of  submodels  to  ade- 
guately  describe  the  engineering  response  to  energetic  loading.  The  far 
field  is  characterized  by  low  strains  and  involves  nominally  elastic  behav¬ 
ior.  Although  this  is  a  simpler  region  to  model,  its  response  is  still  a  func¬ 
tion  of  the  soil  properties  and  reguires  a  material- specific  treatment.  Addi¬ 
tionally,  our  previous  work  (Cole  2015;  Cole  et  al.  2013;  Cole  and  Ketchum 
2013;  Cole  and  Peters  2007,  2008)  clearly  indicates  that  the  location  of  the 
boundary  between  these  two  regions  is  determined  by  the  microscale 
properties  of  the  soil  and  can  be  determined  a  posteriori  from  the  numeri¬ 
cal  solution. 
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The  proper  formulation  of  the  interpartide  contact  behaviors  is  important 
for  accurate  representation  of  the  bulk  properties.  For  example,  the  fric¬ 
tion  coeffident  used  to  calculate  the  tangential  stresses  that  result  from 
the  sliding  action  at  the  particle-  to-  partide  interface  governs  the  transition 
between  elastic  and  plastic  deformation  of  the  soil  skeleton.  The  contad 
equations  that  are  used  in  the  DEM  soil  model  are  described  below.  When 
two  grains  are  found  to  be  in  contad,  a  plane  is  placed  at  the  point  of  con- 
tad  tangent  to  the  surface  of  each  grain.  The  force  between  grains  has  a 
component  normal  to  the  plane  and  a  tangential  component  that  lies  in 
the  plane.  The  normal  and  tangential  components  are  calculated  from  a 
modified  Hertzian  model  with  parallel  viscous  damping  and  a  Coulomb 
fridion  cap  on  the  tangential  force.  The  force  components  depend  on  the 
overlap  between  dilated  grains.  The  overlap  <5  between  a  pair  of  grains  is 
defined  as 


8  =  \d\  —  2R  <  0  (2) 

where  d  is  the  magnitude  of  the  vedor  that  defines  the  distance  between 
the  two  undilated  grains  and  R  is  the  dilating  radius  of  grains  1  and  2.  Di¬ 
lation  of  a  polyhedral  grain  produces  a  similar  polyhedral  grain  with 
curved  edges  and  vertices  with  radii  of  curvature  egual  to  the  dilation  ra¬ 
dius. 

For  the  normal  component  of  the  contad  force,  we  use  a  Hertzian  contad 
force  model  similar  to  the  one  used  by  Lin  and  Ng  ( 1997)  that  was  in  turn 
based  on  the  work  of  J  ohnson  ( 1987) .  The  Hertzian  model  is  based  on  a 
nonlinear  force- displacement  relationship  that  accounts  for  the  increase  in 
contad  area  between  grains  as  the  normal  load  increases.  The  Hertzian 
contad  force  model  was  originally  developed  for  spherical  grains.  The 
eguation  for  the  normal  contad  force,  Fn,  as  a  fundion  of  the  grain  over¬ 
lap,  8,  is 


Fn 


2  G 


3(1  -v) 


-/l1/2s3/2 


(3) 


where 


G  =  the  shear  modulus  of  the  spheres, 
v  =  the  Poisson  ratio  of  the  spheres,  and 
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Re  =  R1R2/  (Ri  +  R2),  where  the  radius  for  each  grain,  Ri  and  R2,  is 
equal  to  the  radius  of  a  sphere  of  volume  equal  to  the  grain 
volume. 


Coulombic  frictional  forces  act  between  each  grain  pair.  The  frictional  con¬ 
tact  force.  Ft,  at  time  m  is  calculated  incrementally  from  the  force  at  the 
previous  time  step  m  -  1  as 


pr  =  771-1 


ktAt  I 


(Vi  -Vi-nj 


(4) 


where 

kt  =  the  tangential  contact  stiffness, 
n  =  the  contact  normal  direction,  and 
At  =  the  time  step. 

The  dependence  of  the  tangential  stiffness  on  the  normal  force  (Lin  and 
Ng  1997)  is 


2[3G\l-v)FnRe]1/3  (5) 

kt  =  - ~ - ■ 

2  —  v 

The  tangential  force  is  also  damped.  Once  the  frictional  force.  Ft,  exceeds 
the  Coulomb  limit,  sliding  at  the  grain  interface  begins.  This  can  be  written 
as  follows: 


if 


p  m+l 

Ft 


>  ixFvm+1 , then 


p  m+l 

Ft 


=  k-Fn 


m+l 


(6) 


where  11  is  the  coefficient  of  friction.  The  tangential  contact  force  (using 
Equation  [4])  is  calculated  incrementally  to  account  for  changes  in  the  di¬ 
rection  of  the  tangential  component  of  the  relative  velocity  and  changes 
that  the  normal  contact  force  produces  in  the  tangential  stiffness. 


After  the  contact  and  body  forces  on  each  soil  grain  are  calculated,  the 
equations  of  motion  are  solved  for  new  positions  and  velocities;  and  time 
advanced  one  step  using  a  Leapfrog  integration  scheme,  which  is  second 
order  accurate  in  time  and  is  conservative.  During  the  simulation,  the  time 
step  is  adjusted  dynamically  using  the  equation 
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n 


(7) 


where  M  is  the  mass  of  the  smaller  of  the  two  grains  in  a  contact  and  the 
effective  contact  stiffness,  keff,  is 


keff  ~ 


2V2G 
3(1  -  v) 


Re1/2S1/2. 


(8) 


During  each  simulation,  the  energy  balance  is  calculated  to  verify  the  self- 
consistency  of  the  simulation.  The  components  of  the  energy  balance  are 
work  performed  on  the  grains,  the  change  in  the  kinetic  and  potential  en¬ 
ergy  of  the  grains,  inelastic  dissipation,  frictional  dissipation,  and  viscous 
drag  losses. 

3.2  Triaxial  shear  experiment 

During  proj  ectile  penetration,  energy  is  transferred  from  the  proj  ectile 
through  the  target  material  via  wave  propagation.  For  unconsolidated 
granular  media,  wave  propagation  is  fundamentally  the  transfer  of  kinetic 
energy  from  one  particle  to  another  through  compression  and  shear.  The 
bulk  material  properties,  such  as  the  elastic  and  shear  moduli,  are  simply 
the  collective  product  of  a  number  of  particle  interactions  and,  therefore, 
can  be  estimated  a  priori  with  the  formulation  of  appropriate  contact  laws 
that  incorporate  the  effects  of  liquid  water  to  account  for  moisture  in  the 
target  material. 

Work  on  the  contact  behavior  of  dry  grains  (Cole  et  al.  2013)  has  shown 
that  the  initial  state  of  force  at  the  contacts  due  to  the  soil-fabric  micro¬ 
structure  and  overburden  pressure  determine  the  dynamic  stiffness  and 
attenuation  for  that  material.  The  near- surface  wave  propagation  and  at¬ 
tenuation  in  moist  soils  depend  not  only  on  these  factors  but  also  on  satu¬ 
ration  level,  pore  geometry,  and  encapsulated  air.  To  examine  the  effects 
of  varying  levels  of  moisture  on  the  dynamic  behavior  of  a  soil,  we  fabri¬ 
cated  a  small-scale  triaxial  shear  test.  The  results  from  these  experiments 
were  used  to  inform  the  DEM  contact  model  development  and  subsequent 
calibration.  As  described  above,  the  forces  at  the  particle-to- particle  inter¬ 
face  in  the  DEM  model  are  parameterized  by  a  shear  modulus,  friction  co¬ 
efficient,  and  damping  coefficient.  We  can  then  use  the  experimental  re- 
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suits  of  the  triaxial  shear  experiments  to  estimate  the  value  of  these  pa¬ 
rameters,  plug  these  values  into  the  numerical  model,  and  use  the  differ¬ 
ence  between  the  simulated  behavior  and  the  observed  behavior  to  im¬ 
prove  our  parameter  estimation.  This  workflow  can  then  be  repeated  for 
varying  soil  conditions  (i.e.,  moisture  contents)  to  parameterize  the  soil's 
mechanical  properties  by  moisture  content. 


3.2.1  Experimental  setup  and  procedure 

The  triaxial  shear  test  setup  was  composed  of  a  linear  piezoelectric  actua¬ 
tor  that  cyclically  loaded  a  granular  sample  and  a  load  cell  that  recorded 
the  sample's  response  to  loading  (Figure  1) . 


Figure  1.  Diagram  of  triaxial  shear  test 
setup  with  a  granular  sample.  This  entire 
setup  sits  inside  of  a  pressure  vessel. 


Actuator 


Pinch 

Collar 

End-Cap 


Granular 

Sample 


End-Cap 
Load  Cell 


Test  samples  consisted  of  American  Society  for  Testing  and  Materials 
(ASTM)  20-30  standard  test  sand  contained  in  latex  membranes,  which 
maintained  the  sample's  form  and  moisture  content  while  allowing  the 
vessel's  pressure  to  confine  the  granular  material.  ASTM  20-30  sand  is  a 
smooth- grained  guartz  sand  that  is  predominantly  graded  to  pass  through 
a  850  pm  (No.  20)  sieve  and  be  retained  by  a  600  pm  (No.  30)  sieve 
(ASTM  International  2013) .  Samples  were  approximately  36  mm  in  diam- 
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eter  and  85  mm  in  length  and  were  moistened  with  a  25%  potassium  io¬ 
dide  solution  (by  weight).  The  ends  of  the  samples  were  enclosed  by  alumi¬ 
num  platens  that  connected  to  the  actuator  and  load  cell.  Rubber  O- rings 
sealed  the  latex  membrane  to  the  aluminum  platens. 

The  actuator's  loading  signal  followed  a  haversine  produced  by  a  com¬ 
puter-controlled  function  generator.  The  loading  signal  was  limited  so  that 
the  resultant  sample  displacement  was  less  than  500  nm  peak-to-peak  and 
the  frequency  was  less  than  150  Hz  (i.e.,  cycles  per  second) .  The  actuator 
contained  a  strain  gauge,  which  measured  the  amount  of  strain  experi¬ 
enced  by  the  granular  sample. 

Sample  moisture  contents  varied  between  dry  and  fully  saturated  (-10%) 
(Note:  for  these  initial  results,  the  "dry"  sample  material  was  not  baked  to 
remove  all  moisture) .  Samples  were  formed  within  a  mold  to  ensure  con¬ 
sistent  sample  shape  and  size  across  all  tests.  Samples  were  prepared  by 
measuring  the  amount  of  material  needed  to  fill  the  mold  and  then  adding 
the  appropriate  amount  of  water  to  achieve  each  desired  moisture  content 
level.  The  material  was  mixed  with  the  water,  inserted  into  the  latex  mem¬ 
brane,  compacted  into  the  mold,  and  then  placed  inside  the  testing  appa¬ 
ratus. 

Once  the  samples  were  properly  aligned  in  the  setup,  the  mold  was  re¬ 
moved;  and  a  standard  preloading  scheme  was  applied  to  initialize  the  ma¬ 
terial  to  a  state  that  was  consistent,  or  as  close  to  consistent  as  possible, 
across  all  samples.  Once  the  standard  initial  state  was  achieved,  the  sam¬ 
ple  was  cyclically  compressed  in  the  axial  direction  for  frequencies  ranging 
from  10  to  150  Hz  for  a  constant  amplitude.  The  actuator  was  controlled 
using  a  displacement- following  feedback  loop  to  ensure  that  the  amplitude 
was  in  fact  held  constant  over  all  the  forcing  frequencies  tested. 

3.2.2  System  dynamics  model 

The  contact  between  DEM  grains  is  modeled  using  a  modified  Hertzian 
contact  model,  which  is  a  nonlinear  viscoelastic  contact  model.  The  system 
dynamics  of  the  triaxial  shear  experimental  setup  can  be  modeled  in  a  sim¬ 
ilar  manner  as  the  DEM  model  to  provide  an  analog  between  the  numeri¬ 
cal  and  physical  experiments.  Assuming  frequency- invariant  stiffness  and 
damping  of  the  experimental  system,  the  measured  force,  F,  is 
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F  —  kx  +  bx  (9) 

where  k,  b,  x,  and  x  denote  the  stiffness,  damping,  displacement,  and  dis¬ 
placement  rate,  respectively.  The  sample  is  subjected  to  a  harmonic  dis¬ 
placement  such  that  the  load  and  displacement  take  the  following  forms: 

F  —  F0cos(cL)t)  (10) 


and 


X  -  X0COS  (cot  +  0) 
x  =  — x0<Dsin(a)t  +  0) 


(11) 


where  0  is  the  phase  lag  and  co  is  the  angular  frequency,  which  is  equal  to 
2 ji  times  the  load  signal  frequency,  f.  Fo  and  xo  are  load  and  displacement 
amplitudes.  A  phase  difference  between  the  displacement  and  load  signals 
is  assumed.  Substituting  Equations  (10)  and  (11)  into  (9)  returns 


F0cos(a)t)  —  /cx0cos(&)t  +  0)  —  bx0(jL>sin(oot  +  0),  (12) 

which  can  be  simplified  with  the  following  trigonometric  identities. 


cos(a  +  /?)  =  cos(a)cos(/?)  —  sin(a)sin(/?) 
sin(a  +  /?)  =  sin(a)cos(/?)  +  cos(a)sin(/?), 
and  rearranged  such  that 


(13) 


[F0  —  kx0  cos (0)  +  bxQod  sin(0)]cos(a)t)  -1-  [/cxosin(0)  (14) 

+  f>xo6jcos(0)]sin(<i>t)  =  0. 

For  Equation  ( 14)  to  hold  true,  the  bracketed  terms  must  both  be  inde¬ 
pendently  equal  to  zero: 


F0  —  /cxocos(0)  +  bx0co  sin(0)  =  0 
/cxosin(0)  +  bx0a)cos(<p )  =  0 


(15) 


These  can  then  be  used  to  solve  for  k  and  b: 
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F0  co s(0) 


Xn 


b  =  — tan(0) 

CO 


(16) 


To  extract  the  elastic  and  viscous  components  of  the  granular  samples,  we 
must  first  determine  the  dynamic  contributions  from  the  experimental 
setup  itself.  The  system's  stiffness  and  damping  were  measured  by  placing 
a  steel  sample  within  the  setup  and  using  Equation  ( 16)  to  calculate  k  and 
b  for  the  measured  results.  It  was  assumed  that  the  steel  sample's  stiffness 
and  damping  were  much  larger  and  lower,  respectively,  than  those  of  the 
system. 


We  then  assumed  that  the  stiffness  and  damping  components  for  the  sys¬ 
tem  combine  linearly  with  those  of  the  sample  such  that  the  total  system 
dynamic  response  is 

Fsys  =  Os  +  km(/))x  +  Os  +  (17) 

where  the  subscripts  s  and  m  refer  to  the  experimental  setup  and  granular 
material  responses  to  dynamic  loading,  respectively.  This  way,  the  dy¬ 
namic  responses  of  both  the  granular  sample  and  experimental  setup  can 
be  separated  from  one  another. 

3.3  Projectile-drop  test 

Ultimately,  the  numerical  modeling  developed  in  this  effort  will  be  used  to 
predict  the  penetration  behavior  and  final  depth  for  a  given  projectile,  im¬ 
pact  conditions,  and  soil  conditions.  A  companion  effort  was  used  to  cali¬ 
brate  the  DEM  model  and  to  validate  the  DEM  modeling  approach,  which 
involved  several  drop  tests  into  dry  and  wet  uncompacted  ASTM  20-30 
test  sand  using  3- dimensional  printed  replicas  of  a  57  mm  M70  projectile. 
Two  projectiles  were  used  for  the  drop  tests:  one  full-sized  projectile  meas¬ 
uring  170  mm  in  length,  56.78  mm  in  diameter,  and  weighing  409.64  g 
and  one  half  sized  projectile  85.05  mm  in  length,  28.35  mm  in  diameter, 
and  weighing  51. 19  g. 


The  impact  material  was  contained  in  a  box  measuring  30  cm  long,  15  cm 
wide,  and  30  cm  high.  The  projectiles  were  dropped  from  a  height  of  24  in. 
for  each  test  and  were  allowed  to  free  fall  into  the  sand  bed.  For  the  dry 
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sand  tests,  the  diameter  of  the  indention  was  measured  and  recorded;  then 
the  impact  crater  was  filled  with  new  dry  sand  and  smoothed  over  between 
each  drop.  For  the  wet  sand  tests,  water  was  added  to  the  drop  area  by  us¬ 
ing  a  bottle  mister  before  each  drop.  The  diameter  of  the  indention  was 
measured  and  recorded;  then  a  sample  of  sand  was  taken  and  weighed  to 
determine  the  moisture  content  at  the  impact  site  for  each  drop.  The  drop 
tests  were  video  recorded  at  389  frames  per  second  so  that  the  proj  ectile 
position  and  velocity  as  a  function  of  time  could  be  estimated  directly  from 
the  frame-by-frame  displacement. 

3.4  Coupling  the  finite  element  method  with  the  discrete  element 
method 

To  reduce  the  overall  computational  complexity  of  the  modeling  effort,  a 
continuum  domain  was  used  to  approximate  the  discrete  particle  domain 
in  the  far  field.  This  is  done  because  the  computational  cost  of  modeling 
the  entire  domain  with  only  discrete  particles  would  be  very  computation¬ 
ally  expensive.  In  the  near  field,  the  soil  was  modeled  using  the  DEM  and 
approximated  as  spherical  particles.  The  size  of  the  domain  that  contains 
the  spherical  particles  was  chosen  to  be  large  enough  such  that  penetration 
of  the  projectile  into  the  soil  is  not  influenced  by  the  continuum  boundary 
(e.g.,  influence  from  reflecting  pressure  waves  at  the  boundary).  To  model 
the  continuum  in  the  far  field,  the  open  source  deal.II  (Bangerth  et  al. 
2007)  finite  element  method  (FEM)  library  is  used.  The  two  types  of  mod¬ 
els  (FEM  and  DEM)  interact  through  a  coupling  surface,  which  is  defined 
on  any  part  of  the  finite  element  domain  that  could  come  into  contact  with 
particles.  The  next  subsections  describe  the  continuum  domain  physics, 
closest  contact  point  determination,  and  interaction  physics  at  the  cou¬ 
pling  surface,  which  are  all  pertinent  to  the  coupling  of  the  finite  element 
and  discrete  element  codes. 

3.4.1  Quasi-static  linear  elasticity 

In  the  far- field  domain  of  the  penetration  model,  we  approximate  the  soil 
as  a  linear  elastic  continuum  and  use  FEM  to  solve  the  system.  Linear  elas¬ 
ticity  is  a  poor  representative  material  model  for  a  real  soil,  but  we  have 
used  a  linear  elastic  constitutive  model  in  lieu  of  a  constitutive  model  that 
includes  plasticity  (e.g.,  Mohr-Coulomb)  to  simplify  the  coupling  scheme 
and  to  reduce  computational  cost  in  this  early  developmental  stage. 
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The  equations  of  motion  that  represent  the  physics  of  the  continuum  do¬ 
main  can  be  seen  in  Equation  ( 18)  where  a  is  the  Cauchy  stress,  F  is  a  body 
force  (e.g.,  gravity),  p  is  density,  and  a  is  the  acceleration.  We  assume  that 
the  influence  of  inertia,  or  dynamics,  in  this  problem  is  small;  so  the  accel¬ 
eration  is  set  to  zero: 


\7-a  +  F  =  pa  =  0  (18) 

We  further  assume  that  in  the  continuum,  the  soil  behaves  as  an  isotropic 
linear  elastic  body  according  to  the  constitutive  law 


o  —  C\s  (19) 

where  C  is  the  fourth  order  elasticity  tensor  and  e  is  the  strain.  We  assume 
that  small- strain  theory  is  appropriate  for  the  continuum  domain  as  we  ex¬ 
pect  the  deformation  to  be  small  in  this  region.  Using  the  small- strain  as¬ 
sumption,  the  strain  tensor  is 


1 ,  ,  (20) 
e  =  -  (Vu  +  (Vu)r) 

where  u  is  the  displacement  and  T  denotes  the  transpose  of  the  gradient 
tensor  Vu.  To  solve  Equation  ( 18),  we  multiply  the  equation  by  a  test  func¬ 
tion,  w,  and  integrate  by  parts  to  get  the  weak  form  of  the  equations  of  mo¬ 
tion 


*  r  r 

Vw:  adVL  —  w  ■  FdD.  + 

J o  J o  Jr 


w  ■  tdF 


(21) 


where  the  integral  over  F  is  the  part  of  the  domain  where  a  traction  t  is  ap¬ 
plied.  For  our  purposes,  the  traction  t  will  be  the  load  applied  by  a  spheri¬ 
cal  particle  when  it  contacts  the  coupling  surface  on  the  continuum  do¬ 
main.  Dirichlet,  or  fixed,  boundary  conditions  are  applied  to  the  outside  of 
the  continuum  domain. 


3.4.2  Closest  contact  point  calculation 

As  particles  travel  through  the  domain,  they  will  interact  with  each  other 
and  they  will  interact  with  the  deformable  continuum  domain  through  the 
coupling  surface.  To  determine  if  a  particle  is  in  contact  with  the  contin¬ 
uum  domain,  the  closest  point  between  the  spherical  particle  center  and 
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the  coupling  surface  needs  to  be  calculated.  Because  our  continuum  do¬ 
main  can  deform  in  time  and  our  finite  elements  at  the  coupling  surface  do 
not  have  to  be  planar,  finding  the  closest  point  poses  a  challenge.  To  ad¬ 
dress  this  problem,  we  solve  the  constrained  optimization  problem 


minz  =  || y(X)  —  P(x)||2  subject  to  X  E  (X1,X2,1)  and  X1,X2  (22) 

e  [-1,1] 


where  z  is  the  norm  sguare  distance  between  a  point  y  on  the  coupling  sur¬ 
face  and  a  particle  center,  P.  Figure  2  shows  the  two  points  and  the  dis¬ 
tance  between  them,  indicated  by  a  dotted  red  line.  When  Eguation  (22)  is 
minimized,  the  result  will  be  a  point  y  on  the  coupling  surface  that  is  the 
closest  point  to  the  spherical  particle. 

Figure  2.  Mapping  between  the  coupling-surface  finite  elements  in  the  real  space 
and  reference  space  during  closest-point  determination.  Point  Pis  the  spherical 
particle  center,  point  ^is  the  closest  point  on  the  reference  element  surface  to  the 
particle  point  P,  and  y(X)  is  the  mapping  of  point  A'onto  the  real  space. 
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To  avoid  the  difficulties  of  finding  the  closest  contact  point  on  the  de¬ 
formed  element,  we  pose  that  there  exists  a  point  on  a  reference  element 
that  can  be  mapped  onto  the  real  space.  The  mapping  from  point  X  in  the 
reference  space  to  point  y  in  the  real  space  is 


(23) 


where 
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Nj  =  the  Lagrange  interpolation  functions, 

Xj  =  the  coordinates  of  the  element  nodes  in  the  real  space,  and 
I  =  the  element  node  number. 

For  the  three  dimensional  problem  we  consider,  the  summation  over  I 
runs  from  1  to  8.  In  our  optimization  problem,  we  constrain  our  minimiza¬ 
tion  statement  by  ensuring  that  the  point  X  in  the  reference  space  lies  in 
the  top  surface  (Xx  <=  [— l,  l],  X2  e  [— l,  l],  X3  =  1).  To  find  the  minimum  of 
Equation  (22),  we  take  the  gradient  of  z  and  set  it  equal  to  zero 


Vxz  =  2 Vxy  ■  (y  -  P)  =  0;  (24) 

and  we  iteratively  find  the  solution,  X,  to  this  equation  using  the  Newton- 
Raphson  method 


xn+l  =xn-  (ylzY^X2  (25) 

where  n  refers  to  an  iteration  step  and  Vx  implies  that  we  are  taking  the 
gradient  with  respect  to  the  reference  coordinates.  After  the  iteration  pro¬ 
cedure  and  the  residual  is  minimized,  if  the  closest  point  is  located  within 
the  bounds  of  the  element,  nothing  has  to  be  done;  and  the  solution  is 
minimized.  If  the  closest  point  lies  outside  the  bounds  of  the  reference  do¬ 
main,  an  additional  procedure  is  performed  that  snaps  the  point  X  to  the 
nearest  point,  which  could  be  on  a  line  segment  or  a  comer.  In  our  code, 
the  closest  point  is  calculated  for  each  element  in  a  subset  of  the  total  cou¬ 
pling  surface  for  computational  efficiency.  The  results  of  such  a  procedure 
can  been  seen  in  Figure  3.  For  this  effort,  we  consider  single  point- wise  in¬ 
teraction  per  particle,  so  the  closest  particle-  surface  pair  is  stored.  If  the 
closest- point  vector  has  a  magnitude  larger  than  the  particle  radius,  noth¬ 
ing  has  to  be  done;  and  the  simulation  marches  forward.  On  the  other 
hand,  if  the  magnitude  of  the  closest- point  vector  is  smaller  than  the  ra¬ 
dius  of  the  particle,  contact  is  detected;  and  an  additional  procedure  needs 
to  be  completed,  which  is  explained  in  the  next  subsection. 
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Figure  3.  Potential  results  of  an  element-wise  closest-point 
calculation  on  the  coupling  surface  with  the  red  segmented  line 
showing  the  closest  point. 


3.4.3  Finite  element  and  discrete  particle  interaction  physics  at  the 
coupling  surface 

The  coupling  surface  is  the  portion  of  the  continuum  domain  that  couples 
the  movement  of  the  spherical  particles  to  the  deformation  of  the  contin¬ 
uum.  At  this  interface,  forces  are  transferred  to  and  from  each  domain 
during  particle-  continuum  contact,  which  is  detected  when  the  magnitude 
of  the  closest- point  vector  is  less  than  the  particle  radius.  As  shown  in  Fig¬ 
ure  4,  the  particle  is  penetrating  the  finite  element  at  a  distance  of  5,  which 
can  be  calculated  by  taking  the  difference  between  the  magnitude  of  the 
closest- point  vector  and  the  particle  radius. 

Figure  4.  A  spherical  particle  penetrating  a 
deformable  continuum  finite  element  at  the 
coupling  surface.  The  nearest  contact  point  is 
indicated  by  the  dashed  red  line,  the  depth  of 
indentation  5  is  indicated  by  the  dashed  blue 
line,  and  /?  is  the  particle  radius. 
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The  contact  force  generated  on  the  continuum  surface  during  impact  is  de¬ 
fined  as 


Pi  =  K8t 


(26) 


where 


Pi  =  the  force  from  the  particle  impact, 

K  =  the  stiffness  of  the  particle, 

&  =  the  penetration  distance,  and 
i  =  the  particle  number  in  the  interaction. 

The  number  of  particles  interacting  with  the  boundary  can  vary  through¬ 
out  a  simulation.  The  particle  impact  force  Pi  is  applied  to  the  continuum 
body  through  the  eguation 


F  = 


w  •  5(x  —  Pi)dT 


(27) 


where 


F  =  the  force  vector  acting  on  the  coupling  surface, 
w  =  the  test  function, 

8(x  -  pi)  =  a  delta  function, 

pi  =  the  closest  contact  point  on  the  coupling- surface  element,  and 

the  sum  is  over  the  number  of  particle  contacts.  This  vector  F  and  the  body 
force  vector  are  added  together  and  are  what  drive  the  deformation  of  the 
continuum  body.  On  the  DEM  side  of  the  code,  an  egual  and  opposite 
force  (to  the  one  applied  to  the  coupling  surface)  is  applied  to  the  particle 
to  drive  the  acceleration  of  the  particle  away  from  the  boundary. 
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4  Results  and  Discussion 

4.1  Triaxial  shear  test 

Recall  that  the  particle- to- particle  interactions  can  be  described  by  the 
grain- scale  elastic  modulus  (or  shear  modulus,  which  is  a  function  of  the 
elastic  modulus  and  Poisson's  ratio),  viscous  damping  coefficient,  and  the 
friction  coefficient.  A  granular  material's  bulk  properties,  such  as  elastic 
modulus  and  viscous  damping  coefficient,  are  a  result  of  the  individual 
particle- to-partide  interactions  in  aggregate.  Therefore,  if  we  can  extract 
the  bulk  elastic  modulus  and  the  damping  coefficient  from  a  sample's  me¬ 
chanical  response  to  small- strain  cyclical  loading  by  using  the  triaxial 
shear  test  and  the  dynamical  system  formulation  described  in  Section 
3.2.2,  then  we  can  then  infer  the  particle- scale  parameters  that  govern  the 
bulk  constitutive  behavior.  We  highlight  that  the  soil  samples  were  tested 
with  no  applied  confining  stress,  resulting  in  weaker  coupling  between 
neighboring  soil  grains. 

In  this  set  of  experiments,  the  forcing  amplitude  (i.e.,  the  maximum  axial 
strain)  was  kept  constant  for  freguendes  ranging  from  10  to  150  Hz.  Fig¬ 
ure  5,  the  triaxial  shear  test  results,  shows  that  the  sample's  spring  stiff¬ 
ness,  which  is  directly  proportional  to  the  elastic  modulus,  generally  de¬ 
creases  as  the  fordng  freguency  increases  (i.e.,  the  material  exhibits 
"softer"  elastic  behavior  where  it  does  not  recover  to  its  original  state  as 
guickly) .  Because  the  particle  velodty  increases  with  freguency,  the  rela¬ 
tive  velodty  between  partides  also  increases  and  leads  to  sliding  at  the  in¬ 
terface.  This  scenario  appears  to  be  corroborated  by  the  decrease  in  stiff¬ 
ness  when  moisture  is  added  to  the  system  where  the  water  is  lubricating 
the  particle  interface,  allowing  more  sliding  to  occur. 

The  viscous  damping  coeffident  also  decreases  with  increasing  freguency 
(Figure  5)  until  reaching  a  nearly  constant  value  for  freguendes  above  SO¬ 
TO  Hz,  which  at  first  is  a  bit  puzzling  if  sliding  is  the  mechanism  for  the 
stiffness  decrease  with  increasing  freguency  (i.e.,  strain  rate) .  But  in  our 
viscoelastic  model  of  the  sample  dynamics  (Eguation  [  16]),  the  damping 
coeffident,  b,  is  proportional  to  the  stiffness,  k.  Therefore,  if  the  ratio  of 
the  tangent  of  the  phase  lag  angle,  tan  (0),  and  the  fordng  freguency,  co,  is 
either  constant  or  slowly  varying  with  freguency,  then  the  damping  coeffi¬ 
dent  would  decrease. 


ERDC/CRREL  TR-17-12 


22 


Figure  5.  Granular  sample  stiffness  (A)  and  damping  coefficients  (b)  for  various 
initial  moisture  contents.  The  material  stiffness  and  damping  coefficients 
decrease  with  increasing  frequency  and  moisture  content.  For  the  frequency 
range  considered,  the  damping  coefficient  does  not  vary  with  moisture  content 
in  a  in  a  discernable  manner;  however,  the  stiffness  transitions  to  a  state  with 
higher  magnitude  when  the  moisture  content  is  negligible. 


4.2  Projectile-drop  test— physical  experiments 

The  results  of  these  projectile- drop  experiments  will  provide  the  data 
needed  to  validate  the  numerical  model  effort.  As  described  previously,  the 
projectile  was  dropped  from  a  height  of  0.61  m  from  the  sand  bed  surface. 
Figures  6  and  7  illustrate  the  projectile  just  before  penetration,  at  midim¬ 
pact,  and  when  the  projectile  has  reached  its  maximum  depth  for  dry  and 
moist  conditions.  The  impact  velocities  (i.e.,  the  projectile  velocities  at  t  = 
0.0  in  Figure  8)  range  from  3. 1  to  3.5  m/  s,  which  is  in  relatively  good 
agreement  with  the  theoretical  impact  velocity  of  3.5  m/  s  of  a  body  accel¬ 
erating  due  to  gravity  from  an  initial  height  of  0.61  m. 

For  each  size  projectile  (full-  and  half-scale  projectile  models)  and  parti- 
de-bed  conditions  (dry  and  2.3%  moisture  content  by  weight),  the  pene¬ 
tration  process  can  be  divided  into  three  stages.  In  the  first  stage,  there  is  a 
roll  off  of  the  velodty  (as  seen  in  Figure  8)  as  the  projectile  first  enters  the 
particle  bed;  and  the  deceleration  is  gradual  as  the  particles  at  the  surface 
provide  little  resistance  to  the  entering  projectile.  As  the  projectile  travels 
farther  into  the  partide  medium,  the  deceleration  reaches  a  maximum  (the 
location  of  the  peak  impad  force  in  Figure  9) .  The  final  stage  is  the  period 


ERDC/CRREL  TR-17-12 


23 


after  the  maximum  impact  force  has  been  reached  until  the  projectile  has 
come  to  rest.  The  magnitude  difference  in  impact  force  between  the  full- 
scale  and  half- scale  projectiles  is  predominantly  due  to  the  mass  difference 
between  the  two  sizes,  where  the  half- scale  projectile  mass  is  13%  the  value 
of  the  full-scale  model. 


Figure  6.  Projectile  model  entering  dry  sand.  The  images 
are  taken  from  the  moment  of  impact  (t=  0),  at  mid¬ 
impact,  and  when  the  projectile  has  reached  its  final  depth. 

(a)  t  =  0  seconds 


(b)  t=  0.2  seconds 


(c)  t=  0.4  seconds 


ERDC/CRREL  TR-17-12 


24 


Figure  7.  Projectile  entering  moist  sand.  Similar  to  the  dry 
sand  images,  the  frames  are  taken  from  the  moment  of 
impact  {t  =  0),  at  mid-impact,  and  when  the  projectile  has 
reached  its  final  depth.  Note  how  the  crater  formation  is 
less  pronounced  as  compared  to  the  dry  particle  bed. 

(a)  t  =  0  seconds 
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Figure  9.  Measured  impact  forces. 


Time  [s] 

There  is  no  observable  difference  in  penetration  behavior  between  the  dry 
and  wet  cases  for  both  the  full-  and  half- scale  models,  likely  due  to  the 
projectile's  low  impact  velocity,  therefore  kinetic  energy,  and  the  small 
change  in  soil  bulk  density  with  the  addition  of  moisture  at  2.3%  by 
weight. 
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4.3  Projectile-drop  test— numerical  model 

The  numerical  companion  to  the  penetration  drop  test  described  above 
was  performed  to  demonstrate  the  feasibility  of  coupling  laboratory  meas¬ 
urements  with  numerical  model  development  and  calibration.  Table  1  lists 
the  simulation  parameters.  We  selected  the  numerical  simulation  values  to 
ensure  stability  rather  than  a  one-to-one  match.  Future  work  will  include 
improvements  to  the  numerical  model  to  enable  direct  comparison  be¬ 
tween  the  model  and  experiments. 

The  preliminary  numerical  model  results  compare  well  in  a  qualitative 
sense  with  the  drop- test  measurements.  The  temporal  evolution  of  the  im¬ 
pact  forces,  plotted  for  the  numerical  model  in  Figure  10,  exhibits  the 
rapid  initial  rise  in  contact  forces  as  the  projectile  enters  the  particle  bed 
and  reaches  a  peak  contact  force  of  350  kN  after  which  the  impact  force 
decreases  at  a  rate  slower  than  the  initial  rise.  The  impact  forces  in  the  nu¬ 
merical  simulation  show  pronounced  high-frequency  fluctuations  that  are 
due  to  the  relatively  small  number  of  grains  that  are  in  contact  with  the 
projectile  at  any  one  time,  as  seen  Figure  11,  which  amplifies  the  effect  that 
any  one  grain  has  on  the  projectile.  We  do  observe  the  development  of  a 
pressure  bulb  around  the  projectile  tip  (in  Figure  11)  although  not  as  pro¬ 
nounced  due  to  the  small  number  of  grains  in  the  simulation.  Again,  we 
anticipate  correlation  between  numerical  and  physical  experiments  to  im¬ 
prove  as  we  increase  the  number  of  grains  while  also  bringing  the  mechan¬ 
ical  properties  of  the  simulated  material  to  match  the  real  material. 


Table  1.  DEM  simulation  parameters. 


Parameter 

Value 

Projectile  mass 

10  kg 

Projectile  diameter 

57  mm 

Grain  mass 

1.0  g 

Grain  diameter 

5  mm 

#  of  grains 

50  x  103 

Shear  modulus 

109  Pa 

Poisson’s  ratio 

0.4 

Time  step 

lO-5  s 

Impact  velocity 

2.0  m/s 
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Figure  10.  Time  trace  of  impact  forces  estimated  from  the  DEM  simulation  of  the 
projectile-drop  test.  The  time  evolution  of  the  penetration  forces  on  the  projectile 
compare  well  in  a  qualitative  sense. 


Time  [s] 


Figure  11.  DEM  simulation  of  the  projectile-drop  experiment.  Example  frame  from 
the  DEM  simulation  showing  the  development  of  the  pressure  bulb  at  the  tip  of  the 
projectile.  The  contact-force-magnitude  units  are  Newtons. 
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5  Conclusions  and  Implications  for  Future 
Research 

The  work  presented  in  this  report  demonstrates  a  framework  that  would 
improve  proj  ectile  penetration  models  by  enabling  the  development  of  bet¬ 
ter  constitutive  models  that  depend  on  only  the  parent  material  character¬ 
istics  (elastic  modulus,  Poisson's  ratio),  the  grain  geometry,  the  friction  co¬ 
efficient,  and  the  volume  of  pore  water.  Of  these  parameters,  the  friction 
coefficient  is  the  least  known  for  a  particular  soil  type.  The  biaxial  shear 
test  experiments,  which  characterize  the  dynamic  behavior  and  bansfer  of 
energy  through  the  soil  fabric,  serve  as  a  good  source  of  data  to  calibrate 
the  friction  coefficient.  Comparison  with  the  drop  test,  and  eventually  gas 
gun  tests,  are  the  basis  for  model  validation. 

5.1  Triaxial  shear  test— numerical  model 

The  frichon  coefficient  calibrahon  would  be  through  direct  comparison  be¬ 
tween  the  physical  experiments  and  numerical  model  simulations  of  the 
same  particle  configurations  and  forcing  conditions.  The  DEM  model  par- 
bde  configurabons  would  be  generated  from  three-dimensional  microCT 
imaging  of  a  soil  sample  of  interest  (Figure  12)  and  subject  to  the  same  un¬ 
steady  forcing  to  calibrate  the  frichon  coeffident  such  that  the  dynamic 
behaviors  match. 

Figure  12.  Using  three-dimensional  microCT  imaging  to  generate  realistic  particle  geometries. 
CRREL  in-house  capabilities  include  a  high-energy  microCT  scanner  that  is  capable  of  imaging 
resolutions  up  to  10  ym.  Segmentation  methods  are  in  development  to  extract  particle 
geometries  as  triangulated  surface  meshes,  as  seen  in  (a),  as  well  as  pore  water  distribution 

and  liquid  bridge  geometries. 

(a)  Detail  of  particle  surface  discretization.  (b)  Sample  of  several  particle  geometries 


generated  from  microCT  scans. 
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5.2  Coupled  DEM  and  FEM  penetration  model 

As  tills  work  continues,  we  plan  to  couple  a  nonlinear,  finite- strain,  ex¬ 
plicit-dynamics,  finite  element  code  to  the  discrete  element  code.  This  will 
give  us  the  ability  to  capture  more  dynamic  and  impulsive  loads  on  a  con¬ 
tinuum  body.  We  also  plan  to  extend  our  constitutive  modeling  capabili¬ 
ties  to  include  nonlinear  material  models  such  as  J  2  plasticity.  As  we  refine 
our  code,  we  continue  to  look  at  ways  to  improve  the  computational  speed 
of  the  coupling- surface  identification,  closest- point  determination,  and 
contact- detection  algorithms  to  be  able  to  model  larger  domains. 

5.3  Final  thoughts 

The  advantage  of  a  micromechanical  approach  to  the  projectile  penetra¬ 
tion  problem  is  the  potential  predictive  power  once  the  material's  parent 
materials  and  interaction  behavior  (e.g.,  moisture  effects  and  friction  be¬ 
tween  grains)  are  determined.  The  projectile  penetration-depth  prediction 
then  moves  beyond  empirical  correlations  that  are  tightly  tied  to  the  con¬ 
ditions  and  soil  characteristics  of  the  test,  thus  lacking  generality,  to  a 
model  that  can  predict  penetration  depth  based  on  guantities  that  are  eas¬ 
ily  measured  in  a  controlled  laboratory  setting.  This  approach  will  ulti¬ 
mately  reduce  the  uncertainty  in  penetration- depth  predictions  and  in 
turn  reduce  MEC  remediation  costs. 
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